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ABSTRACT 

We construct detailed vertical structure models of irradiated accretion disks around 
T Tauri stars with interstellar medium dust uniformly mixed with gas. The dependence 
of the structure and emission properties on mass accretion rate, viscosity parameter, 
and disk radius is explored using these models. The theoretical spectral energy 
distributions (SEDs) and images for all inclinations are compared with observations 
of the entire population of Classical T Tauri stars (CTTS) and Class I objects in 
Taurus. In particular, we find that the median near-infrared fluxes can be explained 
within the errors with the most recent values for the median accretion rates for CTTS. 
We further show that the majority of the Class I sources in Taurus cannot be Class 
II sources viewed edge-on because they are too luminous and their colors would be 
consistent with disks seen only in a narrow range of inclinations. Our models appear 
to be too geometrically thick at large radii, as suggested by: (a) larger far-infrared 
disk emission than in the typical SEDs of T Tauri stars; (b) wider dark dust lanes in 
the model images than in the images of HH30 and HK Tau/c; and (c) larger predicted 
number of stars extincted by edge-on disks than consistent with current surveys. The 
large thickness of the model is a consequence of the assumption that dust and gas are 
well-mixed, suggesting that some degree of dust settling may be required to explain 
the observations. 

Subject headings: Physical data and processes: accretion, accretion disks — stars: 
circumstellar matter, formation, pre-main sequence 



1. Introduction 

With the advent of IRAS, it became clear by the mid-1980s that the infrared excesses of many 
T Tauri stars were probably produced in dusty circumstellar disks rather than gaseous envelopes 
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(e.g., Rucinski 1985). However, it was also recognized that the spectra of the infrared excess of T 
Tauri stars (Rydgren & Vrba 1987) was not of the form predicted either for a steady accretion disk 
(Lynden-Bell & Pringle 1974) or for a flat disk heated by radiation from the central star (Friedjung 
1985), with the implication that the outer disk temperatures are much hotter than predicted by 
the standard models. A variety of mechanisms were invoked to explain this result, but the most 
robust model has proved to be that of the irradiated "flared disk" (Kenyon & Hartmann 1987, 
hereafter KH87). Because the scale height if in a Keplerian (geometrically thin) disk varies with 
the sound speed c s and radius R as H/R oc c s -R 1//2 , any disk temperature distribution which 
decreases less rapidly with radius than i? -1 (expected on quite general grounds) causes the disk 
to become proportionately thicker with increasing R. This "flaring" of the disk allows it to 
absorb more radiation from the star, especially at large distances, resulting in spectral energy 
distributions (SEDs) that reproduce observations much better (KH87; Kenyon & Hartmann 1995, 
hereafter KH95). 

More detailed models of flared disks, which relax the assumption of vertical isothermality 
adopted by KH87, have been constructed by Chiang & Goldreich (1997, hereafter CG), who 
presented a two-isothermal-layer model, to allow for irradiation heating of the upper layers, and 
by D'Alessio et al. (1998, hereafter DCCL), who solve the detailed vertical hydrostatic disk 
structure with irradiation and viscous heating self-consistently. Although the resulting model 
SEDs have been shown to agree with observations of a single T Tauri star in each of these papers 
(not the same star in both), they have not been shown to agree with typical disk emission spectra. 
Moreover, there are reasons to suspect that these models, which assume well-mixed gas and dust 
throughout the disk, may not be able to reproduce many results. In particular, the disk models 
presented in these studies are so vertically thick that, if inclinations were randomly distributed, 
roughly half of all T Tauri stars should be hidden behind their disks for typical radii. Not only 
would this mean large incompleteness in present surveys of star-forming regions, it would also 
suggest that many of the Class I objects, which were originally thought to be protostars with 
infalling dusty envelopes, might simply be T Tauri stars seen edge-on. On the other hand, the 
results of Miyake and Nakagawa (1995) suggested that some T Tauri SEDs could be explained 
with flat disks, possibly implying dust settling to the disk midplane. 

In this paper we construct detailed irradiated disk models for T Tauri stars and compare 
them with typical T Tauri SEDs. This comparison suggests that T Tauri disks are thinner than 
predicted by the models with complete mixing of ISM dust and gas. Survey incompleteness due to 
stars hidden behind their edge-on disks does not seem to be very large. We also argue that most 
Class I objects in Taurus cannot simply be edge-on disks, but are probably protostars surrounded 
by their infalling envelopes, as previously proposed (Adams, Lada, & Shu 1987; Kenyon, Calvet & 
Hartmann 1993). The models suggest that the assumptions of complete mixing of dust and gas 
and standard ISM dust may need to be relaxed, and we explore such new models in a subsequent 
paper. 
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2. 



Disk Models 



2.1. Assumptions and Methods 



We solve self-consistently the complete set of vertical structure equations, including irradiation 
and viscous heating, resulting in detailed profiles of temperature and density with vertical height. 
The model of the disk vertical structure is calculated with the equations and method described 
by DCCL, as summarized briefly here. The disk is assumed to be geometrically thin and in 



The viscosity is described using the a prescription (Shakura & Sunyaev 1972), with a viscosity 
coefficient written as v% = ac s H, where c s is the local sound speed, H is the local scale height of 
the gas and a is the viscosity parameter, assumed constant through the disk. The scale height is 
given by H = c s (T)/fi, where f2 is the angular Keplerian velocity and T is the local temperature. 
A characteristic scale height of the disk is H c , calculated with the sound speed evaluated at the 
midplane temperature T c . The disk has azimuthal symmetry, and its equations are written in 
cylindrical coordinates (R,z), where R is the radial distance from the star, in a direction parallel 
to the disk midplane, and z is the vertical distance from the midplane, in a direction parallel to 
the rotation axis. 

The radiation field is separated into two frequency ranges as proposed by Calvet et al. (1991, 
1992), one characteristic of the stellar effective temperature and the other corresponding to the 
local disk temperature. We use mean opacities to quantify the disk interaction with the two 
different radiation fields. The Planck mean absorption coefficient for a stellar effective temperature 
T* = 4000 K is similar to the monochromatic absorption coefficient at A = 1 pm (e.g. KH87). The 
stellar radiation field is described as a beam of parallel rays impinging on the disk surface at an 
angle cos -1 \i§ relative to the disk normal. 

The differential equations for the disk vertical structure are the same as those in DCCL. 
We write them here in a slightly different and clearer way. First, we consider the equations for 
radiative transfer in the "disk-frequency range" (Mihalas 1978), 



where Fd and Jd are the flux and the mean intensity, respectively, at the disk frequency range, p is 
the mass density, Kp and \R are the Planck mean and the Rosseland mean opacities, respectively, 
and ctr is the Stefan-Boltzmann constant. 

The disk energy balance equation is 



steady state with a constant mass accretion rate M, and in vertical hydrostatic equilibrium. 
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where F is the total flux given by the sum of the radiative, convective, and turbulent energy fluxes 
(DCCL). Here, T v i s is the heating due to viscous dissipation (e.g., Frank et al. 1992), Ti on is 
the heating due to ionization by cosmic rays and radioactive decay (Nakano &: Umebayashi 1986; 
Stepinski 1992) and rV r is the heating by irradiation. This last heating function can be written 
as (see DCCL): 

T irr = 4ttk* p p[J s + Ji] , (4) 

where k* p is a mean opacity calculated using the Planck function at the stellar effective temperature 
T* as a weighting function, and Jj and J s are the mean intensities of the direct and scattered 
stellar radiation field. To calculate the transfer of the stellar radiation through the disk, we use 
the formulation of Calvet et al. (1991, 1992), with the assumption of constant albedo in the region 
of the atmosphere where most of the stellar radiation is scattered. The mean intensity of the 
direct stellar radiation field is given by 

J% — Jirr& ' 5 (5) 

where t sz is the vertical optical depth to the stellar radiation and J{ rr is the mean intensity of the 
stellar radiation at the top of the disk. The scattered mean intensity is 

7 _ s F irr (2 + 3p ) s z _ 3p sF irr z/M . 

Js 47r[l + (2« ? /3)](l-<^2 ) 4vr(l-^) ' W 

where s is the mean albedo, g = [3(1 — s)] 1 ^ 2 , F{ rr is the stellar flux intercepted by the "irradiation 
surface" of the disk. We take po as the cosine of the angle between the incident direction and the 
normal to this surface. The irradiation surface z s (R) is defined by requiring that the optical depth 
in the direction to the star is t s = 2/3, with t s calculated integrating the opacity in the radial 
direction, i.e., 

Xs(R,z) p(R,z) dr , (7) 



where r = (z 2 + i? 2 ) 1 / 2 , and Xs is a mean opacity calculated including true absorption and 
scattering, using the Planck function evaluated at T* as the weighting function. Since this optical 
depth depends on the disk structure, which in turn depends on the irradiation flux, we iterate to 
find the self-consistent disk structure. Notice that in equations (||) and @ we have used the plane 
parallel approximation, so r s ~ t s z /fio, but this approximation is not used to find the height of 
the irradiation surface. 

The turbulent and/or convective energy transport are included as 

— = -V(F rad , F, T,P) -g z p , (8) 

where V is calculated taking into account that the disk may have a finite optical depth (i.e., 
without using the diffusion approximation) and that the convective elements lose energy by 
radiation and by the turbulent flux (see DCCL and D'Alessio (1996) for details). 
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Finally, we have the vertical hydrostatic equilibrium equation 

dP 

where P is the gas pressure, g z is the z component of the stellar gravity, i.e., g z = GM*zj (R 2 +z 2 ) 3 ^ 2 , 
where M* is the stellar mass, and we have neglected radiation pressure and the disk self gravity. 

The boundary conditions for these equations are basically the same as in DCCL. At the disk 
midplane the fluxes are zero because of reflection symmetry. The height of the disk Zoo is an 
unknown boundary where we specify the ambient gas pressure Poo, the emergent fluxes produced 
by viscous dissipation F v i s , ionization heating Fi on (by cosmic rays and radioactive decay) and 
irradiation Fj rr , and the mean radiative intensity at the disk frequency Jd = hFi rr , where h is a 
function of s and fj,o (cf. DCCL). The turbulent and convective energy fluxes equal zero at z^ by 
definition. We assume P^ = 1CT 10 dyn cm -2 , which is an arbitrary value, but small enough that 
its precise value does not significantly affect the resulting disk structure. 

The intercepted irradiation flux Fi rr is calculated as the component of the incoming stellar 
radiant energy normal to the irradiation surface z s (instead of the surface where the pressure is 
the ambient pressure z^, as in DCCL). Note that the normal to the disk surface is not precisely 
in the z direction. Since the slope of this irradiation surface is used to evaluate Fi rr and /io (see 
KH87) we estimate this slope using the following procedure. We fit a power law to the height of 
the irradiation surface, z s {R) ~ R b , in an interval AlogR ~ 0.55 centered at each radius and the 
slope is estimated as dz s /dR ss bz s /R, using the actual value of z s (R). This smoothing procedure 
accelerates the convergence and avoids numerical instabilities, and is justified since irregularities in 
the shape of the disk surface are damped in times shorter than the viscous timescale (Cunningham 
1976; D'Alessio et al. 1999). The irradiation flux is the fraction of the intercepted stellar flux 
absorbed or scattered by the disk. We estimate this fraction as / = (1 — e~ T3 - totat ), where T S)to tai 
is the total optical depth of the disk to the stellar radiation. We calculate the flux as described 
by KH87, but multiply by / to take into account that there can be optically thin regions, which 
cannot absorb all the intercepted stellar radiation. For typical parameters of T Tauri disks, / ~ 1. 

The disk dust opacity is calculated using optical constants from Draine & Lee (1984) and the 
MRN size distribution (Mathis, Rumpl & Nordsieck 1977). For the gaseous component, we use 
the opacities described by Calvet et al. (1991). In this paper, we assume that gas and dust are 
well mixed and thermally coupled. 



2.2. Fiducial Model 

To illustrate the basic results of the detailed disk calculations, we constructed a fiducial disk 
model using typical parameters of CTTS. The central star is assumed to have a radius R* = 2 R Q , 
a mass M* = 0.5 M , and an effective temperature T* = 4000 K typical of K7-M0 T Tauri 
stars in Taurus (KH95). We adopt a mass accretion rate M = 10 -8 M Q yr _1 and a viscosity 
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parameter a = 0.01, as suggested by recent studies (Gullbring et al. 1998; Hartmann et al. 1998). 
Variations of the outer disk radius are easily incorporated in our formalism; the fiducial model 
has Rd = 100 AU. We also include an inner disk hole, consistent with a magnetospheric radius of 
Rhole = 3 R* (Kenyon et al. 1994a; Kenyon, Yi, & Hartmann 1996; Meyer, Calvet, & Hillenbrand 
1997). 

Figure [l] shows characteristic temperatures, heights, surface density, and mass of the fiducial 
disk model as functions of the distance to the central star. The plotted temperatures are 
the temperature at the midplane T c , the irradiation temperature Tj rr , which is the effective 
temperature corresponding to the irradiation flux Fi rr , the surface temperature To, given by 
the temperature at z = z s , the viscous temperature T v i s , defined as urT^ s = F v i s , and the 
photospheric temperature T p h ot , defined as the temperature at the height where the Rosseland 
mean optical depth is tr(z) = 2/3, so it is only plotted for those regions where the disk total 
vertical Rosseland optical depth is larger than this value. 

A comparison of Tp^ot and T v i s indicates that for R ;> 1 AU irradiation is the main heating 
agent of the disk upper layers, while in the inner disk the viscous flux is significant in comparison 
with the irradiation flux. In this region, T c is larger than T p h t because the disk is optically-thick 
to its own radiation, resulting in effective trapping of the viscously-generated flux. The plateau in 
T C (R) ~ 1600 K is a result of dust sublimation, which acts as a kind of thermostat because the 
gas opacity for T <J 2000 K is very low. As the central temperature rises beyond 1600 K, the dust 
opacity which produces the high central temperature by trapping of viscous energy disappears. 
This results in sublimation of dust near the midplane, while the layers near the photosphere retain 
enough dust opacity to produce a nearly fixed central temperature. 

The total optical depth of the fiducial disk model at its own radiation tr becomes ^ 10 for 
R ^ 5 AU. Beyond R ~ 30 AU, it is smaller than 2/3 and T p hot is no longer defined; however, 
the central temperature dependence on radius is similar to that of the regions where the disk is 
optically thick to its own radiation, as long as the optical depth to the stellar radiation is very 
large. To see the reason for this, consider the limit where the cosmic ray and viscous heating can 
be neglected. In this limit, the energy balance equation (eq. |l|) can be written as 

4nHipp(^-^ J d ^j = T irr . (10) 

using F ~ Fd and equation (||). When the disk is very optically thick to the stellar radiation, the 
high optical depth prevents the stellar radiation from reaching directly the disk midplane and the 
irradiation heating Tj rr — > 0, (eq. Q), so the temperature then depends only upon Jd- Similarly, 
near the midplane, Fd ~ (by reflection symmetry) and equation (g) yields Jd ~ Jd(zoo)- So, Jd 
is approximately constant in height at given radius; using the boundary condition, it is given by 



J d ~ hF irr , (11) 
where h is the relevant Eddington factor (cf. § |2.1| ). Thus energy balance yields a central 
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temperature 

o-T 4 ~ vr/iF irr (12) 

and therefore the temperature has the radial dependence of the fourth root of the irradiation flux, 
in this case T c ~ i?~ 3 / 7 . 

The terms on the left-hand side of the energy balance equation ( JlOj ) represent integrals over 
frequency, and in a more complete treatment can be written as 47rp(«pB — KjJd), where kj 
is the mean opacity weighted by the disk mean intensity. We are assuming that Kp = kj, i.e. 
that the disk radiation field has a spectrum comparable to that of the local Planck function B. 
The characteristic temperature of the hotter upper regions where most of the stellar radiation is 
absorbed, and therefore where most of Jd is produced is roughly Ti rr . Figure la demonstrates that 
Tj rr is only slightly higher than T c from the detailed numerical calculations, so it appears to be a 
reasonable approximation to set Kp = kj in our models. 

Other treatments attempt to include the difference in effective wavelength of the radiation 
emitted in different parts of the disk in an approximate way, but these approaches have other 
difficulties. For instance, CG assumed that the hot upper layer which absorbs the central star's 
light radiates into the disk at the maximum temperature, so that Jd can be characterized by 
a temperature ~ T(z oa ), and this is much hotter than the central temperature. However, this 
approximation overestimates the characteristic temperature of Jd because most of the radiation of 
the upper layers is produced at lower temperatures than T{z 00 ), as discussed above. Moreover, 
the CG treatment differs most from ours in the case where the disk is assumed to be optically 
thick to the hot layer radiation but optically thin to its own radiation; but CG ignore the vertical 
stratification of the radiation field, and for optically-thick transfer, the hot layer radiation will be 
degraded into lower-temperature emission as the energy moves toward the disk midplane. 

In the case where the central layers of the disk are optically thick to the radiation impinging 
upon them but optically thin to their own radiation, a typical assumption is to write the energy 
balance equation as 

F irr = 4r p a R T 4 (13) 

namely, optically thick heating is balanced by optically thin cooling. The temperature of the 
assumed isothermal layer obtained from this equation has a much flatter radial dependence than 
iT 3/7 (Canto et al. 1995; Aikawa et al. 1999). We can obtain this equation by taking the integral 
over the vertical optical depth of the energy transport equation (|l|), from which we can write 

F d ~ ±T P CJ R {T) A - ATTTpJd (14) 

where 

<T) 4 = - r Tir'pfdr'p . (15) 

Tp JO 

With Fd ~ Fi rr , this equation becomes 

F irr (l + ATrhrp) « F irr ~ ±T P a R (T) A (16) 
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for Tp « 1, which is similar to equation (13). This shows that the temperature that is obtained 
by using equation (|l^) is a mean temperature over the vertical structure and not the actual 
physical depth-dependent temperature of the interior. This distinction may be important when 
considering disk properties which are very very sensitive to local conditions as for instance, the 
molecular processes which determine the chemical composition of the disk (Aikawa et al. 1997, 
1999). 

Irradiation, enhanced by the disk curvature, also dominates the thermal structure of the disk 
atmosphere. The surface temperature is higher than the photospheric temperature, since the 
stellar incident energy is mostly deposited in the upper layers of the disk; this, in turn, is due 
to the large dust opacity at wavelengths characteristic of the stellar radiation and to the oblique 
entry of this radiation into the disk (Calvet et al. 1991, 1992; Malbet & Bertout 1992; CG). This 
is the temperature inversion found by Calvet et al. (1991), which can produce molecular and 
silicate bands in emission even when the disk is optically thick. The inversion is present even in 
the regions where the disk is optically thin to its own radiation, so T c < Tq (Figure Ilk 



Figure |l]b shows some characteristic heights of the disk as a function of radius. We have 
plotted the irradiation surface height z s , the photospheric height z p hot (where tr = 2/3), the 
maximum height of the disk z^ (where P = P^, see §|2.1|) and the gas scale height H c , calculated 



using the sound speed at the disk midplane. In particular, z s can be taken as proportional to 
H c for R ;> 4 AU, as in the approximation usually found in the literature, with a proportionality 
constant ~ 5. The height z s corresponds to a flared surface, with dz s /dR > z s /R, at least for 
R <^ 340 AU which is the maximum disk radius we have considered. On the other hand, the 
photospheric height has a maximum beyond which it decreases as the disk becomes optically thin 
to its own radiation. However, since the disk is still optically thick to the stellar radiation in the 
radial direction, z s (R) is flared in regions where Zpfo t bends down. The behavior of z^qi and the 
bump shown by H c led Bell et al. (1997) to conclude that the outer disk was in the shadow of the 
inner disk, ignoring the role of the material at the upper disk atmosphere, which is able to absorb 
and scatter stellar radiation effectively. Our results contradict their claim, since a more realistic 
treatment, allowing for the very different wavelengths at which light is absorbed and emitted by 
the disk, results in no shadowed regions. 

The surface density of mass £ is the integral along the vertical direction of the volumetric 
mass density p(R,z). In steady disk models using the a prescription, £ oc M/ut oc M£l/aT (e.g., 
Frank et al. 1992). With T ~ T c oc R~ l l 2 in the outer disk, £ oc R~ 1 , as shown in Figure |c. The 
surface density flattens towards smaller radii, where the midplane temperature increases in the 
optically thick annuli of the disk. The mass of the disk for the assumed viscosity parameter and 
mass accretion rate is ~ 0.016(i?(i/100AU) M for a disk radius in the range 1 <; Rd <^ 300 AU 
(Figure |l|d), which is similar to the typical masses estimated from mm observations of CTTS 
(Beckwith et al. 1990; Osterloh & Beckwith 1995). 
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2.3. Consistency of the plane-parallel approximation 

The incident stellar flux is calculated assuming there is a well-defined surface z s into which 
the radiation enters. Given the irradiation flux, we solve the transfer of stellar radiation, assuming 
a plane parallel geometry. This may not be a good assumption for the outer regions of the disk, 
where the height of the irradiation surface becomes of the order of the radial distance. To test 
this approximation, we have used the two-dimensional disk structure to calculate a posteriori 
the mean intensity of the direct stellar radiation at each depth Ji(R,z) attenuated by the radial 
optical depth, 

J i = — T~ I exp(-r s )dQ , (17) 
7T 4-7T Jn„ 

where the integration is performed over the solid angle subtended by the star, as seen from 
each point of the disk. The stellar disk is divided into bands parallel to the disk midplane and 
symmetrically distributed with respect to the line between the center of the star and a given point 
of the disk. Ignoring the stellar diffuse field, the local heating by stellar radiation can be evaluated 
as 

T irr! i(R,z) = 4ttk* p pJi , (18) 

where the subindex % refers to the direct incident stellar radiation. Thus, the flux corresponding 
to the irradiation input of energy, which has to emerge at the top of the disk is 



Fi(zoa) = / T irri (R,z)dz . (19) 
Jo 

Figure Q shows the flux Fi(z OQ ) and Fi rr given by the stellar flux intercepted by the irradiation 
surface, absorbed by the disk, calculated using the plane-parallel approximation and the iterative 



method described in §2.1. To make a consistent comparison, the flux Fi(z OD ) is smoothed fitting 
a local power law in intervals of AlogR = 0.55 centered at each radius, which is the same kind 
of smoothing procedure we have done in the iterative plane parallel calculation. Figure || shows 
that both fluxes are very similar. Since one of them results from a ray by ray integration of the 
transfer equation of the stellar radiation, and the other was obtained using the plane-parallel 
approximation, we conclude that the latter is a good approximation to describe the effect of the 
stellar irradiation on the disk. 

The plane-parallel approximation breaks down for the transfer of radiation in the "disk- 
frequency" range in the outer disk, where the relevant optical depths can become small. We 
estimate that this does not greatly affect our calculated temperature distributions, which are 
effectively determined by the balance between the irradiation flux heating, which is relatively well 
known as discussed above, and optically-thin cooling, which is insensitive to geometry. A detailed 
test of this approximation requires two-dimensional radiative transfer and is beyond the scope of 
this paper. 
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2.4. The M — a Parameter Space 

We next explore the effect of changing the disk mass accretion rate and the viscosity parameter 
a. Figures || and |I] show the resulting disk structure for M = 10 -9 , 10 -8 , and 10 -7 M© yr _1 , and 
a = 0.001, 0.01, and 0.1, for a fixed central star with M* = 0.5 M , R* = 2 R Q and T* = 4000 K. 
The panels are organized in such a way that the surface density of the disk increases towards the 
bottom and the right of the figure, as M/a increases. 

Figure[|| shows the radial distribution of characteristic temperatures: T c , Tq, T p h ot , and T v i S 



(defined in §2.2). The lower the mass accretion rate, the smaller the contribution of viscous 
dissipation as a heating mechanism of the disk photosphere, which can be seen by comparing T v i s 
and T p h ot . For M <J 10~ 8 M yr _1 , the disk photospheric temperature distribution is dominated 
by stellar irradiation (upper panels). 

The behavior of the central temperature in the inner disk depends upon the relative 
importance of viscous dissipation and irradiation and upon the optical depth of the disk, which in 
turn oc S oc M/a. As the optical depth increases, the ratio T c /T p h ot increases. In all optically 
thick models, the central temperature shows the plateau at ~ 1600if due to dust sublimation, as 
in the fiducial model (§2.2). Since viscous heating at a given radius increases with mass accretion 
rate, the outer boundary of the dust sublimation region moves outwards with higher M. 

For M = 10 M Q yr _1 and a = 0.001, which corresponds to the densest disk model shown 
in Figure |3|, the midplane temperature becomes very high, reaching T c ~ 10000 — 30000 K at 
R < 0.05AU. The dominant opacity source at this high temperature range is the ionization of H, 
He and metals, and the disk probably is subject to thermal instabilities (Kawazoe & Mineshige 
1993; Bell & Lin 1994). 

In all models, irradiation heating enhanced by the flaring of the surface dominates in the outer 
disk. The radius where irradiation begins to dominate (Tj rr > T v i s ) increases as M/a increases. 
In the lowest surface density cases we have calculated, (M/10 _8 )/(a/0.01) = 0.01 and 0.1 (left 
upper corner in Figure §), the temperature distribution flattens out at large radii. In this regions, 
the disk becomes optically thin to both its own radiation and the stellar radiation. In this case, a 
fraction of direct and diffuse stellar radiation is able to penetrate the disk and reach the midplane, 
so T irr (z = 0) ^ (c.f. §2.2). 

Figure || shows the surface temperature inversion discussed in §2.2, in many models over a 
wide range of radii (To > T p h ot for R ^ 0.03 AU). The upper layer temperature is approximately 
the temperature of optically thin dust heated by stellar radiation only geometrically diluted 
(e.g., Calvet et al. 1992; CG; DCCL) and it is independent of M and a, being sensitive to 
dust properties, the effective temperature and luminosity of the central star, and other cooling 
ingredients different from dust (CO, H2O, etc). 

Figure || shows the mass surface density of the different models, with the fiducial S(i?) 
repeated in each panel as a reference. ^From this plot it is clear that £ scales as M/a. 
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As mentioned above, the surface density of a steady a-disk is £ ~ M/aT c , but in the 
outer regions the midplane temperature is almost independent of M and a, since irradiation 
is the dominant heating mechanism, The disk mass, given in Table 1 for three different 
disk radii, depends on the mass surface density of the outer annuli, and scales roughly as 
M d ~ (ik/100 AU)(M/1(T 8 Moyr-^O.Ol/ajM©. 

Regions in the disk are unstable to axisymmetric gravitational perturbations if the Toomre 
parameter 

« = sl < 2 »> 

(Toomre 1964) is less than unity. It has been suggested that disks cannot maintain such regions 
because they will be unstable to non-axisymmetric gravitational perturbations which effectively 
transfer angular momentum, causing rapid accretion (e.g., Pringle 1981). The most unstable 
regions in our models occur at large radii. Table 1 shows the radius Rq beyond which Q < 1. 
The most critical case is the disk model with M = 10 -7 M Q yr _1 and a = 0.001, which is 
gravitationally unstable for R > 14 AU. 



3. Comparison with Observations 

3.1. Calculation of the Model Emission 

Given the disk structure and assuming an inclination angle i between the disk axis and the 
line of sight, we integrate the radiative transfer equation through the disk along rays parallel to 
the line of sight, in a grid of points at the plane of the sky. To calculate the disk thermal emission 
we integrate 

dI therm 



dZ 



-K vP B v e- T »W , (21) 

§ = ■ <») 

where Z is the coordinate along the ray which is zero at the plane of the sky and increases towards 
the observer, t u (Z) is the monochromatic optical depth, B v is the Planck function evaluated at 
the local temperature, k u is the true absorption coefficient, and \v is the total opacity including 
absorption and scattering. 

The contribution of the stellar radiation scattered by the disk is calculated assuming single 
and isotropic scattering, so we integrate the transfer equation 

^jscatt 

-^^ = -a u W(r)B u (n)exp[-T u , rad -T u (Z)} , (23) 

where a v is the scattering coefficient, W(r) is the geometric dilution factor of the stellar radiation 
which reaches radial distance r, B U (T*) is the Planck function evaluated at the stellar effective 
temperature, and T u . ra( i is the optical depth in the radial direction, between the star (assumed 



-12- 



as a point-like source) and the point r in the disk. For each wavelength we calculate T U)Ta( i{R, z) 
in cylindrical coordinates R, z once, and use this array to interpolate during the integration of 
equation (|23|). We calculate the emergent intensity of the disk thermal radiation at 57 wavelengths 
between 0.55 /im and 20 cm, and since the intensity of the stellar radiation peaks around 1 fim, 
we calculate the scattered light contribution only at 13 wavelengths between 0.55 \x and 4.65 /im. 

The intensity emerging at each ray is crossing the plane of the sky is 

I„ = ll herm + I s u catt + r v e~ T - , (24) 

where /* is the stellar intensity, which is added if the ray intersects the stellar disk, and r* is the 
optical depth towards the central star produced by the disk. 

Finally, we calculate the fluxes and images of disks with arbitrary orientations and at different 
wavelengths, convolving with appropriate Gaussian instrumental responses. More details about 
this procedure and useful integration limits can be found in D'Alessio (1996). We take the 
observed SED of V819 Tau which is a weak emission T Tauri Star (WTTS), scaled to 1 L Q , as the 
typical central star for calculating the system SED. 



3.2. SEDs 

3.2.1. The Median Observed SED 

The principal observational constraints we use in this paper come from the SEDs of T Tauri 
stars in the Taurus- Auriga molecular cloud, taken from the compilation of KH95. Taurus may not 
be the most representative region of star formation, but it is nearby, well-studied, and not heavily 
extincted, so that its stellar population is fairly well characterized (e.g., Hartmann et al. 1991; 
Gomez et al. 1993; Briceho et al. 1997, 1998). 

Analysis of the SED of any particular source is complicated by a number of effects. 
One of the principal concerns must be the presence of binary companions, especially infrared 
bright-companions, such as in T Tau (Ghez et al. 1991), which can alter the interpretation of 
the spectrum dramatically. In some cases the near-infrared spectra can be decomposed between 
the objects (e.g., Simon et al. 1992, 1995; Leinert et al. 1993); however, the IRAS fluxes, which 
are crucial to understanding the (cool) outer disk structure, cannot similarly be separated into 
individual components because of poor spatial resolution. Beyond this, there are many other 
factors - variability, uncertainty in extinction corrections, IRAS sensitivity limits, observational 
errors, effects of (generally unknown) inclination which render the fitting of an individual SED 
non- unique. 

For these reasons we have estimated a median TTS SED, with a measure of the typical range 
of SED properties, for comparison with the disk models. To do this we first selected K5-M2 
stars from KH95 to isolate a reasonably restricted range of stellar effective temperature (this also 
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happens to include a large fraction of the known T Tauri stars in Taurus). We eliminated the 
WTTS, to avoid including objects which have little or no infrared excess emission from disks. 
We also eliminated a few objects for which the data is inadequate, or which may have infrared 
excesses dominated by envelope rather than disk emission (e.g., Calvet et al. 1994). Finally, we 
eliminated objects for which the strong optical veiling makes the intrinsic system colors uncertain, 
and therefore have very uncertain extinctions. These criteria resulted in a final list of 39 objects 
for further study, as listed in Table 2. As described in KH95, we adopt the luminosity estimate 
from the reddening-corrected J magnitude, Lj, as the best estimate of the stellar luminosity 
(Table 2). 

We then corrected the individual SEDs for extinction, using the values given in KH95, and 
normalized all the SEDs at 1.6/iin. This normalization was motivated by the expectation (as 
confirmed below) that most of the disk heating in typical TTS is from absorbing light from 
the central star, not by intrinsic accretion energy generation, in which case the disk luminosity 
should scale with the stellar luminosity. For this sample of stars, which have nearly the same 
effective temperature, it is sufficient to normalize in the H photometric band, near the peak of the 
stellar SED. Adopting a near-infrared rather than optical wavelength for normalization also helps 
minimize the effects of errors in extinction corrections. 

The median fluxes and quartile limits resulting from this procedure are given in Table 3, and 
the median fluxes along with the individual observations normalized at 1.6//m are illustrated in 
Figure [| Substantial scatter is observed at each wavelength, but there is a clearly defined trend 
encompassing most objects. The dot-dashed line is the median SED, constructed independently at 
each wavelength. This median is straightforwardly constructed except at A = lOO^m, where 13 of 
the 39 objects are undetected, and at mm wavelengths, where there is very large incompleteness. 
Figure || shows the two estimates of the median, calculated with and without the upper limits 
at 100 /im; the difference at this wavelength between the two estimates is comparable to the 
quartile scatter in the median. In the subsequent comparison between the observed median and 
the models we adopt the 100 /im value which allows for incompleteness (i.e., calculated with upper 
limits). We emphasize that the mm- wave points are meant to be indicative rather than definitive; 
these fluxes are not used as constraints on disk models in this paper because we do not consider 
changing dust opacities (see below). 

This is not the only way of developing a median SED. We have tried normalizing at the J and 
K photometric bands, but the results are essentially the same. A crucial parameter of our models 
is the ratio of the disk brightness to the stellar brightness, and so normalization at the peak of the 
stellar SED rather than at some wavelength dominated by disk emission seems most appropriate. 

In principle, the median flux approach might not yield the correct SED spectral indices as a 
function of wavelength if the underlying SEDs are not smooth. In practice, this does not seem to 
be a problem. The spectral indices derived from our median spectrum agree well with the average 
spectral indices derived by KH95 by straight averaging of spectral slopes without normalization. 
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It might be objected that our method of normalization (basically to the stellar luminosity) 
introduces additional error in cases where the accretion luminosity is not negligible. In such 
cases, the absolute values of L acc and are important. However, it can be seen from Table 2 
that most of the objects selected have quite similar stellar luminosities; half of the sample have 
—0.27 < logL*/LQ < —0.03, so that normalization introduces relatively small shifts in the data for 
many objects. We think that the advantages of our procedure in establishing a well-determined, 
stable benchmark SED to test models against outweigh the disadvantages of giving individual 
differences reduced weight, especially because irradiation heating probably dominates accretion 
energy release in most of the sample. In any event, since we do not include continuum stars which 
are likely to have high accretion rates, neither our sample nor our method are well-suited to an 
exploration of the extreme range of properties among T Tauri disks. 



3.2.2. Fiducial model SED 

Figure |6| shows the median observed SED compared to the SEDs of the fiducial model and of 
a flat disk with the same stellar and disk parameters, assuming an inclination angle of i = 60°. 
Fluxes are presented for the fiducial model truncated at three different outer radii, 30, 100, and 
300 AU. 

It is evident that the fiducial model, viewed at the median inclination i = 60° expected for a 
random distribution of orientations, is roughly consistent with the SEDs of many T Tauri stars 
from near-infrared wavelengths out to around 100 fim. However, the model fluxes are slightly too 
low around A ~ 3/xm, and somewhat too high in the 20 — 100/^m range, when compared with the 
median SED, though the latter problem may be reduced if disk radii are smaller than the typical 
100 AU estimate. 

The largest systematic problem of the fiducial model occurs at mm wavelengths, where the 
model fluxes are too low by a factor of 10 (for A > 1 mm). In contrast with the disk flux at shorter 
wavelengths, which is dominated by optically-thick regions, much of the disk emission in the mm 
range is produced by optically thin regions, and so is strongly dependent upon the dust opacity 
at long wavelengths. The typical practice in modeling mm emission of T Tauri disks is to adjust 
the dust opacity to larger values than resulting from the Draine and Lee (1984) calculations. A 
frequently- used estimate (e.g., Beckwith et al. 1990) would increase the 1 mm opacity by a factor 
of 10 over what we are using here, and such an opacity increase would bring the model fluxes into 
much better agreement with observations. We defer discussion of matching the long-wavelength 
fluxes to the next paper in this series, where we will consider the effect of changing dust opacities. 

The comparison between the fiducial model and the flat disk model shows that flaring is 
relatively unimportant in determining the flux in the 2 < A < 6 /im wavelength range. To estimate 
the required additional heating of the inner region, we made the experiment of increasing the 
irradiation flux by a factor of 2.5 for R < 0.1 AU; this increases the disk photospheric temperature 
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by a factor of 1.2, and gives rise to a spectrum identical to the observed median SED in the range 
1.22 < A < 6 fim. This required large increase in heating suggests that reducing the dust albedo 
will not explain the discrepancy. Alternatively, an increase of M of a factor ~ 3 would produce 
the requiring heating (see § |3.2.3| ). 

The far-infrared flux increases as the outer disk radius increases because the flaring results 
in more light from the central star to be intercepted by the disk. The large far-IR excesses are a 
consequence of the substantial vertical thicknesses of the disk models. The corresponding result 
is that more lines of sight to the central star are strongly extincted by the disk. For the fiducial 
parameters, an outer disk radius of 300 AU results in substantial extinction of the star by the 
disk; at i = 60°, the star would appear to be extincted by approximately Ay ~ 4, as shown in 
Figure ^. Although in principle this extinction would be removed by the reddening correction 
used in constructing the model SEDs, this value is rather large in comparison with the extinctions 
adopted for most Taurus Class II sources (Table 2). 

To illustrate the disk extinction in more detail, we calculated surfaces z(Ay) of constant Ay 
along the line of sight to the star. It is convenient to present these results in terms of a critical 
angle i c , such that cosi c = fj, c = z(Ay)/[z(Ay) 2 + i?^ 2 ] 1 / 2 . Thus, fi c (Ay,Rd) is the cosine of the 
inclination angle at which an observer would view the central star through Ay magnitudes of 
visual extinction for a disk of outer radius Rd- For inclination angles greater than i c = cos -1 fj, c , 
the extinction is larger, and vice versa. 

Figure shows /j, c for Ay = 4, 10 and 30 for the fiducial model. The curves show why the 
models with outer radii < 100 AU exhibited only small extinction toward the central star at 
i = 60° (/x = 0.5), while the 300 AU outer radius model heavily extincted the star. The fiducial 
model has a substantial vertical thickness (cf . Figure ||b) ; the "disk photosphere" that an observer 
might define at visual wavelengths would be quite thick, with a height of nearly half the cylindrical 
radius at 100 AU. 

The large vertical thickness of the fiducial disk model implies that the central star would 
be heavily extincted at a wide range of inclination angles. To illustrate this in detail, Figure [8| 
shows SEDs of our fiducial model as a function of inclination, fixing Rd = 100 AU. As inclination 
increases, the star first becomes heavily reddened. At large inclination, < 0.3, the SED becomes 
separated into two distinct peaks, with the short-wavelength emission resulting from stellar light 
scattered by the disk into the line of sight. 

Large vertical thickness and corresponding large extinction towards the star for a large range 
of inclination angles are also found by CG and Chiang & Goldreich (1998), who assumed that the 
dust and gas are well mixed as we do in this work. However, the disk thickness and the resulting 
infrared emission are larger in the Chiang & Goldreich models than in ours. The difference 
arises for two reasons. First, we include a finite albedo. Second, as discussed in § |2,2j , the disk 
atmosphere in our calculations has a range of temperatures (similarly to a stellar atmosphere), 
and the temperature where most of the emission of the upper atmospheric levels is produced is 
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lower than the single temperature CG assume for these upper layers. 

Conversely, the isothermal disk models of Miyake & Nakagawa (1995) exhibited thinner disk 
structure than found here. However, even with this property, Miyake & Nakagawa argued that the 
SEDs of a few T Tauri stars could be explained only if their disks were geometrically flat, which 
could be attributed to dust settling toward the disk midplane. 

One other property to note is the presence of emission in the 10 \xm silicate band in the 
low inclinations disk models; this band is produced by the temperature inversion in the disk 
atmosphere (Calvet et al. 1991, 1992), which is more conspicuous in the SED of flared disk models 
than in the flat disk. Both, the flat and the flared models have a temperature inversion in the 
vertical direction (Calvet et al. 1991, 1992), but the smaller the angle between the incidence 
direction and the disk normal, the larger the vertical optical depth of the atmospheric layer where 
the largest fraction of stellar energy is deposited and the more important its contribution to the 
disk SED. As inclination increases, a very strong silicate absorption feature develops (as also is 
shown by Chiang & Goldreich 1998, hereafter CG2). 



3.2.3. SED dependence on M and a 



The models presented in §2.4 cover a range of M compatible with observations of veiling in 
CTTS (Valenti, Basri & Johns 1993; Hartigan, Edwards & Ghandour 1995; Gullbring et al. 1998). 
In this section we show how the SEDs depend on M and a to see if different combinations of these 
parameters give a better fit to the median SED than our fiducial model. 



Figure || shows SEDs of irradiated accretion disks with the same parameters as in §2.4, 
namely, M = 10~ 9 , 10~ 8 , and 10 -7 M yr _1 , and a = 0.001, 0.01, and 0.1, at an inclination 
i = 60°. The near-IR emission (~ 2 — 5/xm) is dominated by optically thick regions at R <^ 1 AU 



(see §2.4). The difference between the models with M = 10 -9 and 10 -8 M yr _1 (upper panels 
in Figure ||) is small because irradiation is the dominant heating source at low accretion rates. 
For these models the near IR flux (at i = 60°) is smaller than the median (approximately at the 
level of the lower quartiles). On the other hand, accretion energy is important in enhancing the 
optically thick emission for the M = 10 -7 M yr _1 , and so these models exhibit near IR fluxes 
larger than the median and the upper quartiles. 

We find that an accretion rate of M = 3 x 10 -8 M yr _1 fit the median observed SED 
in the 2-8 A spectral range reasonably well. This is slightly larger than the median accretion 
rate of M ~ 1 x 10~ 8 Moyr^ 1 found by Gullbring et al. (1998). Note that the actual 
median stellar luminosity is L* ~ 10~ 013 Lq, so the required mass accretion rate would be 
M ~ 10°' 34 x 10~ 8 ~ 2 x 10 -8 M yr _1 . The agreement is reasonably satisfactory considering that 
there may be systematic errors in accretion rates at the factor of two level. Moreover, the spread 
between the quartile fluxes is consistent with a random distribution of inclinations and a quartile 
spread of order ~ ±3 in accretion rates consistent with observations (Hartmann et al. 1998). 
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The mid- to far-infrared fluxes (10/um < A < 100 \xm) are much less sensitive to the accretion 
rate because of the importance of irradiation in outer disk regions. In general, the models in the 
grid also show the same excess in this wavelength region relative to the median SED that the 
fiducial model exhibits. The exception is the lowest accretion rate, largest-a model; in this case 
the 10 — 100/im flux is reduced because the very low mass disk becomes optically thin. However, 
these model parameters do not seem to represent a good solution to the excess mid- to far-IR 
fluxes, because the mm-wave fluxes become extremely small, and it would require an enormous 
increase in the long- wavelength opacity to reconcile the model with observations. 

The flux in the mm range tends to scale as the disk mass (e.g., Beckwith et al. 1990), which 
in our grid of fixed radius models scales as £ oc M/a. Figure |9| shows that models with an M/a 
larger by a factor of 10 with respect to the fiducial model (i.e., M = 10 -8 M yr _1 and a = 0.001 
or M = 10~ 7 Mq yr _1 and a = 0.01) can account for the mm emission without changing the 
dust opacity (cf. discussion in previous subsection). However, the masses of these disks are 
Md ~ 0.13 — 0.14 Mq (Table 1), not far from the limit for gravitational stability (Table 1; §3.2.2). 
In addition, the high M of these disks results in near infrared fluxes much higher than the median 
SED. 

At short wavelengths, the extinction of the stellar radiation produced by the disk is shown 
by the difference between the model fluxes and the median. As M/a increases, so does the disk 
density and thus the Ay towards the central star. To illustrate how the extinction depends on 
the disk parameters, Table 4 lists the cosine of the critical inclination angle, fM c , for Ay = 4 and 
Rd = 100 AU. The excess of extinction of the central star radiation (cf. Table 4) and the excess of 
flux in the range 25 — 100 \im are a consequence of the large vertical thickness of the disk models. 



3.2.4- Images 



Figure |10| shows images of the fiducial disk model for inclination angles such that 
cost = fjL = 0, 0.15, 0.3 and 0.4, and a disk radius Rd = 100 AU. The images are calculated at 
A = 1 fim, convolving with a Gaussian point-spread function with a full width half maximum of 
0.1 arcseconds (i.e., 14 AU at a distance d = 140 pc). The contour levels are separated by a factor 
1.58, which is approximately 0.5 magnitudes. For \i < 0.3 the flux from the star is obscured by the 
intervening disk material, and the scattered light intensity relative to the peak intensity is high. 
For n > 0.3 the stellar radiation dominates the image, and the light scattered by the disk is not 
detectable by contrast. 

The image of the edge-on disk (// = 0) shows two elongated reflection nebulosities, almost 
parallel, separated by a dark lane. The apparent thickness of the disk, measured as the distance 
along the polar axis between the center of the elongated nebulae below and above the disk plane, 
is Az app = 0.64" = 89AU, and the aspect ratio of the image, as defined by thickness divided by 
diameter, is ~ 0.5. It is similar in shape to, but larger in aspect ratio than, the disk model shown 
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by Whitney & Hartmann (1992) (cf. their Figure 6), which has an aspect ratio of ~ 0.3. 

There are images of two known edge-on disks in Taurus with which we can compare the 
model predictions. The image of HH30 (Burrows et al. 1996) shows a clearly flared geometry. 
This geometry is not so apparent in the theoretical image because the extinction due to the dust 
in the outer disk is large, even at large heights, resulting in a very thick dark lane which hides the 
curvature of the irradiation surface. The HH30 disk can be detected out to a radius Rd = 250 AU, 
which is larger than our fiducial model; however, increasing the radius of the disk model does not 
change the shape of the image, which still is not as flared as observed (as shown by Whitney & 
Hartmann 1992, cf. their Figure 7). 

Burrows et al. (1996) found that the (isothermal) disk model which produces the best fit 
of the observed HH30 images has a scale height of 15.5 AU at R = 100 AU (similar values are 
found by Wood et al. 1998). We compare this height to the theoretical local scale height (defined 



in §2.1) evaluated at the upper atmospheric layers, since these are the layers where most of the 
scattering is taking place. The scale height determined from fitting the scattered light is factor of 
2 smaller than the gas scale height of the upper atmosphere of our fiducial model, calculated using 
the temperature To . This suggests that the region of the atmosphere where most of the scattering 
is happening is colder than the upper atmospheric layers of the fiducial model, supporting the idea 
of a less flared irradiation surface. 

HK Tau/c (Stapelfeldt et al. 1998; Koresko 1998) has a shape closer to the image of the 
edge-on model, i.e., two parallel reflection nebulae separated by a dark lane. The apparent 
thickness of this disk is Az app ~ 29 ± 3 AU (Koresko 1998) and the thickness of the dark lane is 
0.13 times the apparent width of the nebulae (Stapelfeldt et al. 1998). In the case of the edge-on 
model, the relative thickness of the dark lane is ~ 3 times the observed one. Again, we conclude 
that the disk model is too thick to explain the observations. The disk model that Stapelfeldt et 
al. (1998) proposed as the best fit for HK Tau/c has a scale height of 3.8 AU at R = 50 AU. The 
scale height of the disk model evaluated at T rr is a factor of 2 higher, and it is even larger at To, 
suggesting again that the disk has to be thinner and colder than the model. 

In summary, the thickness of the predicted dark lane in our edge-on fiducial model is larger 
than what is observed in HK Tau/c, and the overall disk structure is inconsistent with the 
observations of HH30, which also probably require a thinner absorption layer, suggesting that our 
fiducial model is a factor of ~ 2 too thick in the vertical direction. 



4. Discussion 
4.1. Edge-on Disks and Class I sources 



If T Tauri disks have substantial geometrical thicknesses, then a sizeable number of objects 
should be viewed edge-on, through the disk. Many heavily-extincted young stellar objects could 
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be edge-on disk systems, detected only as infrared sources. These considerations may suggest that 
some of the Class I infrared sources, usually interpreted as protostars surrounded by infalling 
envelopes (e.g., Adams et al. 1987; Butner et al. 1991; Kenyon, Calvet & Hartmann 1993) could 
be instead simply T Tauri stars viewed through their disks (CG2). For example, CG2 presented a 
SED calculation for a disk system viewed at an intermediate angle which agrees reasonably well 
with that of an individual Class I source in Taurus. 

As in the case of T Tauri stars, in considering Class I sources it is helpful to analyze the 
properties of an entire population of objects rather than simply fitting an individual source. One 
clear prediction is that edge-on systems should appear less luminous than they really are, since 
the disk emission is reduced at large inclinations and the central star becomes heavily extincted 
(Figure |8|). If Class I objects were simply Class II (face-on disk) systems viewed edge-on, then 
the Class I sources should appear to have systematically lower luminosities than Class II objects. 
Yet the luminosity distributions of Class I and Class II sources in Taurus (taken from KH95) 
are nearly identical (Figure |ll]). To make a quantitative test, we have constructed a predicted 
luminosity distribution for edge-on disk systems as follows. For a given stellar luminosity, L*, we 
calculate a frequency distribution of observed luminosities for the star+disk system assuming a 
random distribution of inclinations. We scaled the fluxes of the fiducial model to L* to obtain the 
luminosity of the star+disk system at a given inclination. We restrict the range of inclinations 
to < \i < 0.3, which corresponds to SEDs with \F\ increasing with A, as expected for Class 
I objects. Finally, we convolved this frequency distribution with the observed distribution of 
stellar luminosities of Class II sources in Taurus. The luminosity distribution resulting from this 



procedure, shown in Figure 11, indicates that Class I sources are brighter than expected from 



edge-on disks, suggesting that only a modest fraction of these sources are likely to be edge-on disk 



systems. The luminosity distribution resulting from this procedure, shown in Figure 11, differs 
substantially from that of the currently recognized Class I sources in Taurus, suggesting that only 
a modest fraction of these sources are likely to be edge-on disk systems. 

Another problem with interpreting Class I sources as edge-on disks is that the observed SEDs 
do not agree with the models over a significant range of inclinations. In particular, our models 
tend to predict extremely large silicate absorption, which is not observed in general (Kenyon, 
Calvet & Hartmann 1993). This can be illustrated by the color-color diagram K — [12] vs. 
[12] — [25], shown in Figure |l~2|. The circles represent the observed colors of Class II (open circles) 
and Class I (dark circles) sources in Taurus (KH95), the line corresponds to the colors of the 
fiducial model for fj, < 0.35. The predicted color K — [12], which quantifies the slope of the SED 
between the near and mid IR, increases with the inclination angle and for [i <^ 0.3 becomes larger 
than the maximum observed K — [12] for a Class I source. Only disks within a narrow range of 
inclination angles would have the IR colors and bolometric luminosities consistent with observed 
Class I sources. Because the probability of observing an object in a given range of inclination is 
proportional to the range in fj,, the model predicts a much larger number of sources with K — [12] 
and [12] — [25] colors much redder than observed. 
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In summary, while our models suggest, in agreement with CG2, that some Class I sources 
may be edge-on disk systems, most Class I sources are unlikely to be edge-on disk systems because 
they are not systematically underluminous in comparison with Class II (face-on disk) systems, and 
because their infrared colors and SEDs are inconsistent with disk model predictions for all but a 
narrow range of inclinations. 



4.2. "Missing" sources 



The model calculations shown in Figure 12 indicate that approximately 30% of all T Tauri 



disk systems should have colors and infrared SEDs which are not characteristic of known objects. 
In Taurus, with a total Class II T Tauri population of ~ 100 objects (KH95), this would correspond 
to about ~ 40 systems with peculiar SEDs. While this appears to be strongly inconsistent with 
observations, one must also consider selection effects. In particular, since edge-on systems are 
fainter, one must consider whether surveys have systematically missed these sources. 

The first question we address is whether edge-on disk systems would be detectable as IRAS 
sources. As above, we scaled the edge-on disk models to the stellar luminosity distribution of Class 
II sources, but consider the flux distribution at 60 fim instead of the luminosity distribution. We 
find that approximately 80 % of the edge-on disks should be detectable, using a flux limit of 0.5 Jy 
at 60 microns as adopted in the survey of faint Taurus sources by Kenyon et al. (1994b). Scaling 
from the known Class II sources in Taurus, this would suggest that there should be approximately 
34 edge-on sources in the region which should be detected from IRAS surveys. 

To test this hypothesis, we use the results of Kenyon et al. (1994b), who searched the IRAS 
Serendipitous Survey Catalog for previously unidentified point sources with fluxes exceeding 0.5 
Jy at 60 microns, adding some other candidates with fluxes above 0.3 Jy at 25 or 60 microns 
provided they were detected in at least two bands from the SSC and the PSC2 catalogs. We 
selected objects from Kenyon et al. (1994b) which are not galaxies as found by near-IR imaging 
and which also were not classified as Class I or II sources. There are 19 such sources, out of which 
we estimated 8 candidates by comparison with theoretical SEDs at high inclinations. Figure [D| 
shows two of these candidates for edge-on systems: 104267+2221 and I04413+2608A. However, 
this estimate is likely to be a generous upper limit to the number of possible edge-on sources, since 
there is a lack of multiwavelength data for most of these objects. Some of these sources could still 
be galaxies (9) as mentioned by Kenyon et al. (1994b). In addition, not all objects with near-IR 
and far-IR peaks are necessarily edge-on disks. It is known that many young binary systems have 
greatly differing extinctions to the individual stars, and this multicomponent SED could mimic 
that of an edge-on disk. But even counting all these candidates and including the ~ 10 Class I 



objects that overlap the predicted distribution of edge-on systems (Figure 11), we end up with 19 
candidates, a factor of ~ 2 lower than expected. 



In summary, our disk models with well-mixed gas and dust appear to predict a larger number 
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of low-luminosity, but detectable, highly reddened objects than consistent with current surveys of 
Taurus. 



5. Conclusions 

We have presented models of irradiated accretion disks with well-mixed ISM dust and 
compared them with observations representative of the entire population of Classical T Tauri stars 
in Taurus. We find that our models can explain the near-infrared fluxes of the median SED of T 
Tauri stars with a disk accretion rate consistent with the mean value estimated by Gullbring et 
al. (1997). However, our models appear to be too geometrically thick at large radii, as suggested 
by: (a) larger far-infrared disk emission than the typical SEDs of T Tauri stars; (b) scattered 
light images with wider dark dust lanes than observed in two objects in Taurus; and (c) model 
predictions of large numbers of T Tauri stars hidden by their edge-on disks which contradict 
current surveys. 

Solar nebula theory predicts that dust settling and coagulation could occur on timescales that 
are comparable to or shorter than T Tauri lifetimes, (e.g., Hayashi, Nakazawa & Nakagawa 1985; 
Weidenschilling Sz Cuzzi 1993; Weidenschilling 1997). Because dust dominates the disk opacity, 
settling can strongly reduce the geometrical thickness of the disk, and may be required to explain 
the observations of T Tauri stars (Miyake & Nakagawa 1995). We will consider this possibility in 
a subsequent paper. 

This work was supported in part by the Visiting Scientist Program of the Smithsonian 
Institution, by NASA grant NAG5-4282, by Academia Mexicana de Ciencia and Fundacion 
Mexico-EEUU para la Ciencia, by DGAPA/UNAM and CONACyT, Mexico. 
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TABLE 1 

MODEL DISK MASSES AND TOOMRE RADIUS 
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TABLE 2 
STARS USED FOR MEDIAN SED 



Object 


Ay 


log L j /L(7) 


HO Tau 


1.15 


-0.72 


HN Tau 


0.53 


-0.66 


DM Tau 


0.00 


-0.60 


GO Tau 


1.21 


-0.55 


FM Tau 


0.71 


-0.49 


FS Tau 


1.89 


-0.49 


DP Tau 


1.49 


-0.39 


IP Tau 


0.25 


-0.37 


CY Tau 


0.09 


-0.33 


V710 Tau 


0.90 


-0.27 


DI Tau 


0.78 


-0.21 


DS Tau 


0.31 


-0.19 


IQ Tau 


1.27 


-0.19 


CoKu Tau/3 


3.35 


-0.18 


DH Tau 


1.27 


-0.17 


HK Tau 


2.39 


-0.17 


DD Tau 


0.78 


-0.14 


DQ Tau 


0.99 


-0.14 


AA Tau 


0.50 


-0.13 


JH 112 


3.32 


-0.13 


LkCa 15 


0.65 


-0.13 


V819 Tau 


1.40 


-0.09 


GM Aur 


0.15 


-0.08 


GI Tau 


0.90 


-0.07 


CI Tau 


1.83 


-0.06 


FY Tau 


3.56 


-0.05 


DN Tau 


0.50 


-0.04 


V955 Tau 


2.76 


-0.03 


BP Tau 


0.50 


-0.02 


GK Tau 


90 


07 


DO Tau 


2.70 


0.08 


Haro 6-37 


2.17 


0.11 


LkHa 332/G1 


3.25 


0.15 


GG Tau 


0.78 


0.18 


DF Tau 


0.22 


0.20 


UZ Tau E 


1.52 


0.20 


UY Aur 


1.40 


0.30 


FV Tau 


4.84 


0.34 


V807 Tau 


0.00 


0.61 
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TABLE 3 
MEDIAN SED AND QUARTILES 



Wavelength (/xm) 


log \F\ (med) 


upper 


lower 


0.36 


-10.01 


-10.24 


-9.70 


0.44 


-9.70 


-9.88 


-9.56 


0.55 


-9.48 


-9.63 


-9.34 


0.64 


-9.28 


-9.43 


-9.17 


0.79 


-9.10 


-9.23 


-9.04 


1.22 


-9.01 


-9.04 


-8.96 


1.63 


-9.00 


-9.00 


-9.00 


2.19 


-9.10 


-9.17 


-9.08 


3 45 


-9 43 


-9 49 


-9 32 


4.75 


-9.60 


-9.72 


-9.50 


10.60 


-9.84 


-10.07 


-9.53 


12.00 


-9.81 


-9.96 


-9.50 


25.00 


-9.89 


-10.10 


-9.71 


60.00 


-10.16 


-10.56 


-9.97 


100.00 


-10.43 


-10.73 


-10.13 


800.00 


-11.81 


-11.89 


-11.58 


1100.00 


-12.26 


-12.33 


-11.86 


1300.00 


-12.59 


-12.81 


-12.30 



Note. - The 100/mi entry is the median flux for the sample including upper limits; the lower 
quartile flux cannot be directly determined because of the large number of non-detections, so the 
quartile spread is assumed to be the same as when considering only detections (see text). 
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TABLE 4 
MODEL DISK CRITICAL ANGLE 



7 71 /r 

logM 


7 

loga 


Hc(A v = 4) 


i c 


-9 


-1 


0.25 


75.5 


-9 


-2 


0.34 


70.1 


-9 


-3 


0.41 


65.8 


-8 


-1 


0.34 


70.1 


-8 


-2 


0.41 


65.8 


-8 


-3 


0.46 


62.6 


-7 


-1 


0.41 


65.8 


-7 


-2 


0.46 


62.6 


-7 


-3 


0.50 


60.0 
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Fig. 1. — Structure of the fiducial model. Panel (a)- characteristic temperatures: T c (solid), 
To = T(z = z s ) (dotted), Tj rr (short-dashes), T v i s (dots and long dashes) and T p h ot (crosses). Panel 
(b)- Characteristic heights: Zoo/R (dotted), z s /R (solid), H c /R (dashed) and z p hot/R (dot-dashed). 
Panel (c) - Mass surface density. Panel (d) - Cumulative mass of the disk as a function of distance 
to the central star. 
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Fig. 2. — Comparison between the irradiation flux from the plane-parallel calculation (solid line) 
and Fi(zoo) given by a ray by ray integration of the transfer equation of the stellar radiation (dotted 
line) . 
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Fig. 3. — Temperatures of disk models with log(M/M yr _1 ) = —9, —8, —7, log a = —3, —2, — 1, 
around a central star with M* = 0.5 M Q , = 2 R© and T* = 4000 K. The plotted temperatures 
are: T c (solid line), To (dotted line), T p hot (crosses) and T v i s (dot-dashed line). 
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Fig. 4. — Mass surface density of the disk models plotted in Figure ||. The fiducial model is repeated 
in each panel (dotted line) as a reference. 
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Fig. 5. — Observed fluxes of Taurus-Auriga pre-main sequence sources (KH95) normalized at 
A = 1.6 fim and median SED. The curves are the median with upper limits at 100 /im (dashed 
line) and the median without upper limits (dot-dashed line). 
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Fig. 6. — SED of the fiducial model at an inclination angle i = 60° relative to the line of sight, for 
three disk radii R d = 30 AU (dot-dashed line), 100 AU (solid line), and 300 AU (dashed line). The 
median observed SED (points) and quartiles (error bars) and the disk model irradiated as a flat 
disk (dotted line) are also shown. 
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Fig. 7. — Cosine of critical inclination angle for the fiducial model, and Ay = 4 (solid line), Ay = 10 
(dashed line) and Ay = 30 (dotted line), as a function of disk radius. 
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Fig. 8. — SEDs of the fiducial disk model with Rd = 100 AU and different inclination angles, 
cos(i) = 0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.75 and 1.0 (from bottom to top). As a 
reference, the SED corresponding to cos(z) = 0.5° is shown with dotted line. The median observed 
SED is also plotted (points and error bars). 
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Fig. 9. — SEDs of the same models shown in Figure ||, for i = 60° and = 100 AU, compared 
with the median of the observations and quartiles (points and error bars, respectively). The 
SED of the fiducial model is repeated in each panel (dotted line) as a reference. In the case of 
M = 10 -7 Mq yr _1 and a = 0.001, the high flux at short-wavelengths is due to the emission of the 
inner wall at the hole radius Rhole = 3i2*. This inner region has a very high temperature as can 
be seen in Figure |3], which probably corresponds to a thermally unstable solution. The SED of a 
disk model with the same parameters but a Rhole = 4i2* does not show the high optical-near IR 
emission. 
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Fig. 10. — Images of the fiducial disk model for different inclination angles i = 90° (// = 0), i = 75.5° 
(fi = 0.25) and i = 60° (fx = 0.5). The maximum contour level is I ma x = 11-56 mJY/beam, 
the minimum is I m i n = 0.03 mJY/beam, and each contour level is a factor 1.58 larger than the 
previous one (this corresponds to variations of 0.5 mag). The dotted contour corresponds to the 
half maximum brightness of each image. The images are calculated assuming a distance of d = 140 
pc to the Taurus molecular cloud. 
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Fig. 11. — Luminosity distributions of Class I and II Sources in Taurus (KH95), and edge-on disk 
models (see text) 
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Fig. 12. — Color-color diagram for K-[12] and [12]-[25]. Observed colors are indicated by open 
circles (Class II Sources) and filled circles (Class I Sources), and are taken from KH95. The solid 
line represents the colors of the fiducial model for inclination angles % > 69.5° (fi < 0.35). 
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Fig. 13. — Examples of candidates for edge-on disks from the sample of Kenyon et al. (1994b). 
SEDs of the fiducial disk model for different inclination angles are shown for reference, open circles 
correspond to upper limits. 



